
options(scipen=20)

library(MASS)
library(WhatIf)
library(RColorBrewer)
library(simcf)
library(verification)
library(foreign)
library(nnet)
library(xtable)
library(apsrtable)
library(ggplot2)
library(dplyr)
library(stargazer)
library(haven)

data <- read_dta("cas_chimet.dta")

### Basic analysis (Models found in Table 1 of manuscript) ###

model1<-reg~custodial+questioned+org+black+latino+other+polint+eff2+educ+fem+young+old+dem+ind+inc1+married+unemp
m1<-glm(reg~custodial+questioned+org+black+latino+other+polint+eff2+educ+fem+young+old+dem+ind+inc1+married+unemp, 
        family="binomial",data=data,weights=hweight)

data1<-data[ which(data$reg==1), ]
model2<-vote~custodial+questioned+org+black+latino+other+polint+eff2+educ+fem+young+old+dem+ind+inc1+married+unemp
m2<-glm(vote~custodial+questioned+org+black+latino+other+polint+eff2+educ+fem+young+old+dem+ind+inc1+married+unemp, 
        family="binomial",data=data1,weights=hweight)


model3<-act_index~custodial+questioned+org+black+latino+other+polint+eff2+educ+fem+young+old+dem+ind+inc1+married+unemp
m3<-glm(act_index~custodial+questioned+org+black+latino+other+polint+eff2+educ+fem+young+old+dem+ind+inc1+married+unemp, 
        family="poisson",data=data,weights=hweight)

library(AER)
dispersiontest(m3)

stargazer(m1,m2,m3, no.space=TRUE,type="text",
          single.row = TRUE,
          covariate.labels=c("Correctional Control","Detained by Police",
                             "CSO Connection","Black","Latino","Other Race",
                             "Political Interest","Political Efficacy","Education",
                             "Female","Age: 18-34","Age: 65+",
                             "Democrat","Independent","Income","Married",
                             "Unemployed","Constant"),
          dep.var.labels.include = FALSE,
          model.numbers=FALSE,
          star.char=c("*","**","***"),
          star.cutoffs = c(0.05,0.01,0.001),
          column.labels = c("Registered to Vote","Voted in 2012","Non-Electoral Participation"),
          dep.var.labels   = "")

### Breaking out the participation battery (coefficients found in Figure 1 of manuscript, tables A6-A7 of Appendix)

m1<-glm(petition~custodial+questioned+org+black+latino+other+polint+eff2+educ+fem+young+old+dem+ind+inc1+married+unemp, 
        family="binomial",data=data1,weights=hweight)

m2<-glm(share~custodial+questioned+org+black+latino+other+polint+eff2+educ+fem+young+old+dem+ind+inc1+married+unemp, 
        family="binomial",data=data1,weights=hweight)

m3<-glm(protest~custodial+questioned+org+black+latino+other+polint+eff2+educ+fem+young+old+dem+ind+inc1+married+unemp, 
        family="binomial",data=data1,weights=hweight)

m4<-glm(letter~custodial+questioned+org+black+latino+other+polint+eff2+educ+fem+young+old+dem+ind+inc1+married+unemp, 
        family="binomial",data=data1,weights=hweight)

m5<-glm(donate~custodial+questioned+org+black+latino+other+polint+eff2+educ+fem+young+old+dem+ind+inc1+married+unemp, 
        family="binomial",data=data1,weights=hweight)

m6<-glm(volunteer~custodial+questioned+org+black+latino+other+polint+eff2+educ+fem+young+old+dem+ind+inc1+married+unemp, 
        family="binomial",data=data1,weights=hweight)

m7<-glm(opinion~custodial+questioned+org+black+latino+other+polint+eff2+educ+fem+young+old+dem+ind+inc1+married+unemp, 
        family="binomial",data=data1,weights=hweight)


stargazer(m1,m2,m3,m4, no.space=TRUE,type="text",
          single.row = TRUE,
          covariate.labels=c("Correctional Control","Detained by Police",
                             "CSO Connection","Black","Latino","Other Race",
                             "Political Interest","Political Efficacy","Education",
                             "Female","Age: 18-34","Age: 65+",
                             "Democrat","Independent","Income","Married",
                             "Unemployed","Constant"),
          dep.var.labels.include = FALSE,
          model.numbers=FALSE,
          star.char=c("*","**","***"),
          star.cutoffs = c(0.05,0.01,0.001),
          column.labels = c("Petition","Share","Protest","Letter"),
          dep.var.labels   = "")


stargazer(m5,m6,m7, no.space=TRUE,type="text",
          single.row = TRUE,
          covariate.labels=c("Correctional Control","Detained by Police",
                             "CSO Connection","Black","Latino","Other Race",
                             "Political Interest","Political Efficacy","Education",
                             "Female","Age: 18-34","Age: 65+",
                             "Democrat","Independent","Income","Married",
                             "Unemployed","Constant"),
          dep.var.labels.include = FALSE,
          model.numbers=FALSE,
          star.char=c("*","**","***"),
          star.cutoffs = c(0.05,0.01,0.001),
          column.labels = c("Donate","Volunteer","Opinion"),
          dep.var.labels   = "")

# plot (Figure 1)

coefs_m1 = as.data.frame(summary(m1)$coefficients[2:4,1:2])
names(coefs_m1)[2]="se"
coefs_m1$vars = c("01","02","03")
coefs_m1$group = "01"

coefs_m2 = as.data.frame(summary(m2)$coefficients[2:4,1:2])
names(coefs_m2)[2]="se"
coefs_m2$vars = c("01","02","03")
coefs_m2$group = "02"

coefs_m3 = as.data.frame(summary(m3)$coefficients[2:4,1:2])
names(coefs_m3)[2]="se"
coefs_m3$vars = c("01","02","03")
coefs_m3$group = "03"

coefs_m4 = as.data.frame(summary(m4)$coefficients[2:4,1:2])
names(coefs_m4)[2]="se"
coefs_m4$vars = c("01","02","03")
coefs_m4$group = "04"

coefs_m5 = as.data.frame(summary(m5)$coefficients[2:4,1:2])
names(coefs_m5)[2]="se"
coefs_m5$vars = c("01","02","03")
coefs_m5$group = "05"

coefs_m6 = as.data.frame(summary(m6)$coefficients[2:4,1:2])
names(coefs_m6)[2]="se"
coefs_m6$vars = c("01","02","03")
coefs_m6$group = "06"

coefs_m7 = as.data.frame(summary(m7)$coefficients[2:4,1:2])
names(coefs_m7)[2]="se"
coefs_m7$vars = c("01","02","03")
coefs_m7$group = "07"


coefs<-rbind(coefs_m1,coefs_m2,coefs_m3,coefs_m4,coefs_m5,
             coefs_m6,coefs_m7)

lab <- c(
  "01" = "Petition",
  "02" = "Share",
  "03"= "Protest",
  "04"= "Letter",
  "05"= "Donate",
  "06"= "Volunteer",
  "07"= "Opinion"
)

p1<-ggplot(coefs, aes(x=vars,y=Estimate, ymin=(Estimate - 1.96*se), 
                     ymax = (Estimate + 1.96*se)))+
  geom_errorbar(data=coefs, mapping=aes(x=vars, ymin=(Estimate - 1.96*se),ymax = (Estimate + 1.96*se)),
                                  position=position_dodge(width=.5),width=0.2, size=1)+ 
  facet_wrap(~group,ncol=7)+facet_grid(. ~ group,labeller=labeller(group=lab))+
  geom_point(data=coefs, mapping=aes(x=vars, y=Estimate), 
             size=2, shape=21, fill="white",position=position_dodge(width=.5))+
  theme_bw(base_size=14,base_family="Times")+geom_hline(yintercept = 0,linetype="dotted")+ 
  theme(strip.background = element_rect(fill="white"))+
  scale_x_discrete(breaks=c("01","02","03"),labels=c("Correctional","Detained","CSO"))+
  theme(strip.text.x = element_text(size = 12))+
  theme(axis.text.x = element_text(angle = 45, hjust = 1))+
  labs(y="Marginal Effect\n",x=" ",title="")

ggsave(p1, filename = paste("coeff_act_items", ".pdf", sep = ""),width = 8, height = 4)


# Moderating effect of CSO connections (Table 2 of manuscript)

m1<-glm(reg~custodial*org+questioned*org+black+latino+other+polint+eff2+educ+fem+young+old+dem+ind+inc1+married+unemp, 
        family="binomial",data=data,weights=hweight)

data1<-data[ which(data$reg==1), ]
m2<-glm(vote~custodial*org+questioned*org+black+latino+other+polint+eff2+educ+fem+young+old+dem+ind+inc1+married+unemp, 
        family="binomial",data=data1,weights=hweight)

m3<-glm(act_index~custodial*org+questioned*org+black+latino+other+polint+eff2+educ+fem+young+old+dem+ind+inc1+married+unemp, 
        family="poisson",data=data,weights=hweight)

stargazer(m1,m2,m3, no.space=TRUE,type="text",
          single.row = TRUE,
          covariate.labels=c("Correctional Control","CSO Connection",
                             "Detained by Police","Black","Latino","Other Race",
                             "Political Interest","Political Efficacy","Education",
                             "Female","Age: 18-34","Age: 65+",
                             "Democrat","Independent","Income","Married",
                             "Unemployed","Correct*CSO","Detained*CSO",
                             "Constant"),
          dep.var.labels.include = FALSE,
          model.numbers=FALSE,
          star.char=c("*","**","***"),
          star.cutoffs = c(0.05,0.01,0.001),
          column.labels = c("Registered to Vote","Voted in 2012","Non-Electoral Participation"),
          dep.var.labels   = "")

# Visualize (Figure 2 of manuscript)
model<-act_index~custodial*org+questioned*org+black+latino+other+polint+eff2+educ+fem+young+old+dem+ind+inc1+married+unemp

sims <- 10000
pe<-c(m3$coefficients, m3$zeta)
vc<-vcov(m3)
simbetas <- mvrnorm(sims,pe,vc)

xhyp <- seq(0,1,1)  
nscen <- length(xhyp)  

none <- questioned <- custodial <- cfMake(model, data, nscen)  
for (i in 1:nscen) {
  none <- cfChange(none, "org", x = xhyp[i], scen = i)
  none <- cfChange(none, "questioned", x = 0, scen = i)
  none <- cfChange(none, "custodial", x = 0, scen = i)
  
  questioned <- cfChange(questioned, "org", x = xhyp[i], scen = i)
  questioned <- cfChange(questioned, "questioned", x = 1, scen = i)
  questioned <- cfChange(questioned, "custodial", x = 0, scen = i)
  
  custodial <- cfChange(custodial, "org", x = xhyp[i], scen = i)
  custodial <- cfChange(custodial, "custodial", x = 1, scen = i)
  #questioned <- cfChange(questioned, "custodial", x = 0, scen = i)
}

y1 <- loglinsimev(none, simbetas, ci=.95) 
y2 <- loglinsimev(questioned, simbetas,ci=.95) 
y3 <- loglinsimev(custodial, simbetas,ci=.95)

pe<-rbind(y1$pe[1],y1$pe[2],
          y2$pe[1],y2$pe[2],
          y3$pe[1],y3$pe[2])

low<-rbind(y1$lower[1],y1$lower[2],
           y2$lower[1],y2$lower[2],
           y3$lower[1],y3$lower[2])

hi<-rbind(y1$upper[1],y1$upper[2],
          y2$upper[1],y2$upper[2],
          y3$upper[1],y3$upper[2])


y<-cbind(pe,low,hi)
y<-data.frame(y)
y<-y %>% mutate(org=c("1","2","1",
                      "2","1","2"))
y<-y %>% mutate(contact=c("0","0","1","1","2","2"))

colnames(y)<-c("pe", "lower", "upper","org","x")

df<-data.frame(y)
grouplab <- c(
  "1" = "No CSO Connection",
  "2" = "CSO Connection")

p<-ggplot(df, aes(x=x,y=pe,group=org, ymin=lower, ymax=upper))
p1<-p+ylim(.5,3)+ geom_point(data=df, mapping=aes(x=x, y=pe,shape=org), size=4,
             position=position_dodge(width=.5))+
  scale_shape_manual(values=c(16,18),labels = c("No CSO Connection","CSO Connection"))+
  geom_errorbar(data=df, mapping=aes(x=x, ymin=upper, ymax=lower),
                position=position_dodge(width=.5),
                width=0.1, size=.4) +
  theme_bw(base_size=16,base_family="Times")+
  theme(panel.border=element_rect(fill=NA, colour=NA), 
        legend.title = element_blank(),
        legend.position="bottom")+
  # panel.grid.minor=element_line(colour=NA))+
  scale_x_discrete(breaks=c(0,1,2),labels=c("None","Detained","Correctional Control"))+
  theme(strip.background = element_rect(fill="azure"))+
  theme(strip.text.x = element_text(size = 10))+
  labs(y="Expected Value of \nNonvoting Participation\n",x=" ",
       title="  ")

ggsave(p1, filename = paste("cas_interaction", ".pdf", sep = ""),width = 7.5, height = 5)






